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Abstract 

Rubella is a completely immunizing and mild infec- 
tion in children. Understanding its behavior is of 
considerable applied importance because of Congen- 
ital Rubella Syndrome, which results from infection 
with rubella during early pregnancy and may entail a 
variety of birth defects. The dynamics of rubella are 
relatively poorly resolved, and appear to show con- 
siderable diversity globally. Here, we investigate the 
behavior of a stochastic seasonally forced susceptible- 
infected-recovered model to characterize the deter- 
minants of these dynamics and illustrate patterns 
by comparison with measles. We perform a sys- 
tematic analysis of spectra of stochastic fluctuations 
around stable attractors of the corresponding deter- 
ministic model and compare them with spectra from 
full stochastic simulations in large populations. This 
approach allows us to quantify the effects of demo- 
graphic stochasticity and to give a coherent picture 
of measles and rubella dynamics, explaining essential 
differences in the patterns exhibited by these diseases. 
We discuss the implications of our findings in the con- 



text of vaccination and changing birth rates as well 
as the persistence of these two childhood infections. 



1 Introduction 

Rubella is a completely immunizing, directly trans- 
mitted infection, generally expressed as a mild, and 
potentially even asymptomatic childhood disease ■ 
As a result, rubella tends to be under-reported, and 
its dynamics are fairly poorly characterized. Never- 
theless, since infection during early pregnancy may 
cause spontaneous abortion or Congenital Rubella 
Syndrome (CRS), which may entail a variety of birth 
defects [l^l , understanding the dynamics of rubella is 
of considerable applied importance. Dynamical fea- 
tures of rubella may alter the CRS burden via their 
effects on the average age of infection. Episodic dy- 
namics may increase the average age of infection, as 
the intervals between larger outbreaks provide the op- 
portunity for individuals to age into later age classes 
j2Cl |53| . Likewise, local extinction dynamics can al- 
low individuals to remain susceptible as they age into 
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their childbearing years [37|, |36[ , resulting in the po- 
tential for a considerable CRS burden once rubella is 
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Figure 1: (Color online) Time series of the case reports of rubella and the corresponding spectrum in (a), 
(d) the Hidalgo district, Mexico; (b), (e) Japan and (c), (f) the province Ontario, Canada. To resolve low- 
frequency periodicities these time series include short intervals of vaccination (years 1998-2000 for Mexico 
and 1989-1992 for Japan). Before the spectrum was taken, each series was normalized, setting the mean to 
zero and the variance to unity. The smooth spectrum (thick green lines) was obtained from the raw spectrum 
(thin gray lines) using two passes of a 3-point moving average of the spectral ordinates. The dashed black 
lines are 90% confidence limits on the smooth spectrum. The method of computation of the spectra and 
confidence limits is described in detail in Chapter 4 of Ref. [31 ■ 
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re- introduced. 

Empirically, rubella seems to be linked to either 
i) annual dynamics, as in Mexico [s^, Peru [s^, or 
arts of Africa [2^; ii) spiky dynamics, as in Canada 
and iii) some hint at multi-annual regularity, as in 
Japan 51 1 , England and Wales Q , and various Euro- 
pean countries [13], see typical time series in Figure 
1. Spectral analysis of time series is particularly use- 
ful for understanding temporal patterns exhibited by 
different data [4l|, [3l ■ The characteristic feature of 
rubella spectra are an annual peak at 1 year and a 
multi-annual peak at 5-6 years exhibited by all data 
in Figure 1. Rubella also seems to experience regu- 
lar fade-outs [s^, which is of great epidemiological 
importance, particularly in the context of increasing 
global control efforts. The propensity for stochastic 
extinction is characterised by the critical community 
size (CCS), below which the infection tends to go 
extinct in epidemic troughs. Analyses of dynamics 
in Mexico and Peru suggest a CCS of over 10^ for 
i5|, |37[ , as compared to the estimated value 
10^ — 3 X 10^ for measles in England and 



rubella 
of 2.5 X 
Wales 0] 

The two key ingredients underlying models of 
childhood diseases like rubella and measles are i) sea- 
sonality in transmission due to schooling patterns and 
ii) demographic stochasticity due to the discrete na- 
ture of population [13, [H, [H H^]. Although var- 
ious approaches have been used to understand the 
dynamics of rubella 0, H, S El, most of the anal- 
ysis has been essentially deterministic. Keeling et. 
al. [131 considered a term-time forced susceptible- 
infected-recovered (SIR) model and compared its dy- 
namics to rubella data in Copenhagen. From this, 
they concluded that the dynamics of rubella may 
result from switching between two cyclic attractors 
(annual and multi-annual limit cycles) of the deter- 
ministic model. Although the deterministic analysis 
they present is comprehensive, there is only a lim- 
ited amount of evidence to suggest that the switch- 
ing will occur in contexts that include demographic 
stochasticity. In particular, in this work [soj stochas- 
ticity was introduced into the model as multiplicative 
noise of arbitrary amplitude instead of using, for in- 
stance, standard stochastic simulations based on the 
Gillespie algorithm Such simulations produce 



the exact dynamics of the system subject to demo- 
graphic stochasticity, which for large populations is 
described by the deterministic model with additive 
noise 



Bauch et. al. [6[ studied a term-time forced 
susceptible-exposed-infected-recovered (SEIR) model 
and showed that frequencies obtained from the linear 
perturbation analysis of the deterministic model are 
in good agreement with positions of the peaks in spec- 
tra of data records of various childhood infections. 
The application of this approach to rubella data for 
Canada predicted two distinct peaks at 1 and 5.1 
years, close to what we see in Figure 1 (f). To date, 
there has been little work overall on a fully stochastic 
approach dealing with demographic stochasticity. 



Here, we use this approach to characterize different 
rubella dynamics and illustrate patterns by compar- 
ison with measles. To this end, we perform the the- 
oretical analysis of spectra of stochastic fluctuations 
around stable attractors of a seasonally forced deter- 
ministic SIR model and compare them with spectra 
obtained from full stochastic simulations based on a 
modification ^ of the algorithm by Gillespie [isj . 
The mathematical techniques employed in this work 
have been developed for ecological and epidemiologi- 
cal models 34, 44, 4^ and applied to model temporal 
patterns of measles and pertussis [i^, [^ [13] . The pic- 
ture that emerges to explain rubella dynamics is close 
to that proposed in |6| but goes beyond it because 
our analysis allows us to obtain the full structure of a 
spectrum (as opposed to the deterministic analysis of 
[6| where only frequencies of the spectral peaks could 
be predicted). By introducing key spectral character- 
istics we systematically investigate how the dominant 
period, amplitude and coherence of stochastic fluctu- 
ations change across a broad range of epidemiological 
parameters. We then discuss the implications of our 
analysis in the context of changing birth rates and 
vaccination levels, as well as their implications for 
the persistence of measles and rubella. 



3 



2 Methods 

2.1 Model 

The individual-based stochastic model we explore in 
this paper follows a standard seasonally forced SIR 
structure 2^, 3|- At any time t, it consists of a dis- 
crete population of constant size N divided into com- 
partments of susceptible, S{t), infected, I{t), and 
recovered, -R(t), individuals. Susceptible individu- 
als become infected (and infectious) at a frequency 
dependent rate f3{t)I{t)/N, where f3{t) is a season- 
ally varying transmission rate. For childhood dis- 
eases, (3{t) captures an increase in the number of 
contacts between school children during terms with 
respect to holidays |47|. We consider a sinusoidally 
forced f3{t) = (3o{l + ecos27rt), where /3o is the aver- 
age transmission rate and e is the amplitude of sea- 
sonal forcing. Infected individuals recover at con- 
stant rate u {1/v is the average infectious period). 
As is common in the mathematical epidemiology lit- 
erature [2^, 0] , we restrict our attention to the case 
when birth and death rates /i is the average life- 
time) are equal, and thus the total population size 
is constant. This allows us to reduce the number 
of independent variables to two and define the state 
of the system as ct = {S{t), I{t)}. From /Sq, v and fi 
we can express one of the most important epidemio- 
logical parameters 0, 01 , the basic reproductive ra- 
tio i?o = /^o/(i^ + A')- ^0 is the average number of 
secondary cases caused by one infectious introduced 
into a fully susceptible population; i?o will be used 
throughout the text. 

2.2 Theoretical Analysis 

Two main approaches can be used to investigate the 
dynamics of the stochastic model formulated above. 
An analytical approach starts from the formulation 
of the model as a master equation for the probability 
of finding the system in state a with S(t) susceptibles 
and I{t) infectives at time t 0, 0, Hj] ■ Much un- 
derstanding about the stochastic dynamics relevant 
for recurrent epidemics can be gained if this equation 
is expanded in powers of \/y/N [H^l. An extensive 
discussion of this approach has already been given 



at length in the literature in the context of epidemic 
models, and we refer the reader to [H H, Q for for- 
mal details. Here we describe only the aspects which 
are important for this paper. In essence, the method 
involves the substitutions S{t) = Ns + y/Nxs{t) and 
I{t) — Ni + y/Nxiit) in the master equation which 
can be then expanded to obtain two systems of equa- 
tions [13] ■ At the leading order, the expansion gives 
rise to a set of ordinary differential equations for 
the mean variables, i.e. the fractions (densities) of 
susceptible and infected individuals, s and i. These 
equations are the same as the standard deterministic 
SIR model with sinusoidal forcing (see e.g. [l^). At 
ncxt-to-leading order it gives rise to a set of stochastic 
differential equations for fluctuations of susceptible, 
xs{t), and infected, xi{t), individuals about the mean 
behavior given by the deterministic model [H^ ] . From 
these equations we are able to analytically calculate 
power spectra of fluctuations for susceptibles, Ps{f), 
and infectives, P/(/), as functions of frequency /. We 
are interested in the endemic behavior of the model 
so the spectra correspond to the fluctuations about 
stable attractors of the deterministic model which for 
e = and e > are the endemic fixed point [l| and 
stable limit cycles with a period that is an integer 
multiple of a year li^l , respectively. Further techni- 
cal details relating to analytical calculations are given 
in the Supplementary Electronic Material. 

In our analysis wc will focus on a spectrum Pi{f), 
which for simplicity will be denoted as P{f), and its 
three characteristics, namely dominant period^ am- 
plification and coherence P, [H^l, see Figure 2. We 
define the dominant period of stochastic fluctuations 
as inverse frequency of the maximum of the highest 
stochastic peak. We also compute the total spectral 
power which equals the area under a power spectrum 
curve. This quantity defines the ability of the system 
to sustain oscillations of all frequencies and shall be 
referred to as the amplification of stochastic fluctu- 
ations. Finally, the coherence is defined as the ratio 
of spectral power lying within 10% from the domi- 
nant period and the total spectral power. It serves 
to measure how well-structured oscillations about the 
dominant period are. 

As a rule, a theoretical spectrum of unforced epi- 
demic models (e = 0) has one peak [l|, |4^, while 
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Figure 2: Schematic plot of a power spectrum of 
stochastic fluctuations for infectives, P{f). The 
quantities used in the comparative analysis of dif- 
ferent spectra are the dominant period (the inverse 
of the main peak's frequency), amplification (total 
area under the power spectrum curve) and coherence 
(ratio of the shaded area to the total area). 

it can have several peaks of different amplitudes for 
e > 0, IH] ■ Away from bifurcation points of the de- 
terministic model, one of them is usually much higher 
than the others. We are not aware of any work as- 
sessing the relevance of secondary peaks to recurrent 
epidemic behavior seen in real data. The highest 
peak, however, has been shown to be important in 
understanding the inter-epidemic periods observed in 
time series of pertussis and measles [13, IS] , and is 
therefore used in the definition of a spectrum's char- 
acteristics in this paper. 

2.3 Simulations 

We simulate the model using an extension of Gille- 
spie's algorithm [3,[25| which produces stochastic tra- 
jectories for {S{t), I{t)} in continuous time. These 
are processed further to compute numerical spec- 
tra and test them against the theoretical prediction 
for P{f). The simulation length is 500 years and 
the first 50 years are discarded. In numerical work 



a time series for fluctuations xi{t) is obtained as 
xi{t) = [I{t) — Ni] /Vn, where i is the fraction of in- 
fectives averaged over many realizations of the model. 
From Xi{t) we compute a spectrum P(/) using the 
discrete Fourier transform. For e > we also present 
a spectrum of the entire 'signal' (scaled by popula- 
tion size A''), I{t), which will be referred to as a full 
spectrum. By definition P(/) includes only stochas- 
tic peaks, while the full spectrum includes both de- 
terministic peaks corresponding to a limit cycle and 
stochastic peaks corresponding to fluctuations about 
it. For each set of parameters, 250 simulations start- 
ing from random initial conditions are recorded, and 
all final spectra are averaged over those where no ex- 
tinctions occurred during 500 years. 

3 Results 

We compare the numerical and theoretical predic- 
tions for different spectra and the three measures we 
have conventionally chosen to characterise them. In 
the beginning we explore a large region of parameter 
space and later discuss the main findings for rubella 
and measles. 

3.1 No seasonal forcing: e = 

3.1.1 Theoretical and simulation results 

We first restrict our attention to the case when there 
is no seasonality, for which an explicit expression for 
the analytical spectrum can be found (see the Sup- 
plementary Electronic Material). The deterministic 
SIR model has only one endemic fixed point provided 
i?o > 1 m, I3l- Figure 3 shows analytical and nu- 
merical results for the dominant period of stochastic 
fluctuations about it as well as their amplification and 
coherence for a range of basic reproductive ratios, Rq, 
and infectious periods, l/v. Analytical spectra can 
be obtained for any i?o > 1- In practice long numer- 
ical simulations may not be feasible for all parame- 
ter combinations where i?o > 1 because the system 
experiences frequent extinctions when the infectious 
period is short. The results presented in Figure 3 are 
for a typical population size of a large city, A^ = 10^, 
and the domain where this happens is shown in gray. 
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Figure 3: (Color online) Analytical (top) and simulation (bottom) results for the (a), (d) dominant period of 
stochastic fluctuations about the endemic fixed point, (b), (e) amplification and (c), (f) coherence. The black 
stairstep graph bounds the gray region where all simulations went extinct within 500 years. Approximate 
parameter values for measles and rubella are depicted as dashed circles. Parameters: e = 0, /i = 0.02 f/y 
and N = f 0^. 
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In the region of parameter space amenable to the ex- 
ploration of the long-term dynamics of the model, 
we observe that the structure of the spectra uncov- 
ered by the theory is clearly visible in the simula- 
tions too. Small systematic deviations between the 
two predictions are expected and occur close to the 
extinction boundary. These are mainly reflected in 
the broadening of a spectrum and appearance of sec- 
ondary peaks. As a consequence, amplification, see 
Figure 3 (e) , is slightly increased (coherence. Figure 3 
(f), is correspondingly decreased) in the simulations 
in the area adjacent to the gray region. In addition to 
these changes, an increase of the dominant period of 
stochastic fluctuations in the simulated spectra may 
be observed [H^, 44 1- This effect will be discussed 



in more detail in the next section where seasonality 
is included, as without seasonality, it is only barely 
apparent (compare Figure 3 (a) and (d)). 

3.1.2 Implications for rubella and measles 
dynamics 

The structure discovered in Figure 3 allows us to de- 
rive an initial picture of the dynamics of rubella rel- 
ative to measles. Based on estimates of parameters 
typical of these diseases for the pre-vaccination era 
0, SO, H, S, B we superimpose their approximate 



locations in all panels of Figure 3. The rubella esti- 
mates are for Mexico and Canada, and the measles 
estimates are for England and Wales. For rubella 
the infectious period, l/v, is about 18 days and Rq 
ranges from 3.4 to 9.5 in Mexico, in particular it was 
estimated to be 5.3 to 6.7 in the Mexico's district 
Hidalgo (Figure 1 (a)). Rq in the Canadian province 
Ontario ranges from 4.6 to 6.5 where the upper bound 
is the estimate for the years shown in Figure 1 (c). 
For measles we have taken the most frequently used 
values for large cities (e.g. London) in England and 
Wales before vaccination, 1/v about 2 weeks and i?o 
around 18. 

The results so far ignore the seasonality of trans- 
mission rate and so are insufficient to explain the 
patterns of measles in which it plays a pivotal role 
28j . However, they have important impli- 



of stochastic fluctuations for this disease is close in 
form to that obtained for the unforced model which 
correctly predicts a dominant period of about 5-6 
years (Figure 3 (a) and (d)) as seen in the data (Fig- 
ure 1). This period is similar to the natural period 
of small amplitude perturbations from the endemic 
fixed point recovered in the purely deterministic set- 
ting [SOj . 

Our analysis of the stochastic model allows us to 
quantify other features of fluctuations using amplifi- 
cation and coherence, see Figure 3 (b)-(c) and (e)-(f). 
For rubella, the amplification is large indicating that 
the epidemic patterns of the unforced model repre- 
sent high amplitude oscillations. High coherence sug- 
gests that only a few of the frequencies involved in the 
stochastic fluctuations account for most of the vari- 
ance of time series. This peculiar type of dynamics 
sets rubella close to the extinction boundary. Large 
coherent multi-annual epidemics with troughs deeper 
than in the region with higher i?o cause regular fade 
outs. In the next section, we discuss how these de- 
scriptions of the stochastic dynamics of rubella are 
changed in the presence of seasonality and compare 
it with the dynamics of measles. 

3.2 Seasonal forcing: e > 

3.2.1 Theoretical and simulation results 

For e > the spectra are associated with stochastic 
fluctuations about stable attractors of the determin- 
istic model, i.e. stable limit cycles of a period in 
multiples of a year [i^, 0] . The seasonally forced de- 
terministic SIR model has a complex bifurcation di- 
agram with reg imes where multiple limit cycles may 
coexist [3l|, |52|. For high birth rates and high sea- 



sonality regions corresponding to chaotic dynamics 



cations for understanding the dynamics of rubella. 
As we shall confirm shortly, for e > the spectrum 



are found 2l|. Across most of the range of param- 
eter space we explore here, either annual or biennial 
limit cycles are present. We performed the theoretical 
analysis for these attractors for different parameters 
and found that the agreement between the theory 
and simulations is in general excellent. Nevertheless, 
small discrepancies are again expected if a limit cy- 
cle is not sufficiently stable and/or a population is 
small. In particular, this happens near the extinc- 
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Figure 4: (Color online) Analytical spectra of stochastic fluctuations around an annual cycle (red lines) and 
full numerical spectra (black and green lines) for the seasonally forced model. The red line is indistinguishable 
from the green line except for in the vicinity of 1 year. For clarity in the comparison one of the peaks at 1 
year is cropped in the main panel while the entire spectra can be seen in the inset. Parameters: e ~ 0.05, 
H = 0.02 1/y, i?o = 4, 1/:/ = 18 d, TV = 10^ (black) and N = 10^ (green). 



tion boundary, and is therefore relevant for rubella. 

For parameters reflecting rubella, an annual limit 
cycle is found in the deterministic model. To illus- 
trate the effects of population size on simulated spec- 
tra we show in Figure 4 an analytical spectrum about 
this attractor (red line) and full numerical spectra for 
N = 10^ (black hnc) and iV = 10^ ( green line). For 
the larger population size, the simulated spectrum 
exhibits two types of peaks. There is a dominant an- 
nual peak (corresponding to the deterministic annual 
cycle) and a subdominant broad multi-annual peak 
(corresponding to stochastic fluctuations about it). 
The latter is indistinguishable from the theoretical 
spectrum (red and green lines coincide except for in 
the vicinity of 1 year). 

The amplitude of the deterministic peak scales 
with population size and for ~ 10^ the peak at 
1 year becomes subdominant, see the black line in 
Figure 4. This indicates that the contribution of an 
annual component in the time scries decreases with 
decreasing N. As for fluctuations beyond the annual 
component, at least two stochastic peaks at 5.8 and 
2.9 years can be clearly seen. Although the theory 
does not capture them in full, the agreement is still 



good and, more importantly, the systematic quali- 
tative changes can be predicted. For small popula- 
tions, the dominant period of fluctuations in simula- 
tions is slightly increased and their variance is dis- 
tributed over a larger range of frequencies with re- 
spect to theoretical predictions 14J]. This is com- 
patible with a general observation of the increased 
stochasticity and therefore much more irregular dy- 
namics in small populations (27l . [sij . The example 
we presented here was for Rq = 4, which neighbors 
the extinction boundary. Simulations for larger Rq 
show smaller deviations from analytical calculations 
even for populations as small as = 10^. 

Another situation which the analytical theory can- 
not fully account for is stochastic switching be- 
tween different attractors of the deterministic model 
[itI . [sot ■ The computation of an analytical spectrum 
about a limit cycle requires the knowledge of its ge- 
ometric orbit (see Supplementary Electronic Mate- 
rial). In the sinusoidally forced SIR model several 
stable attractors coexist in the regions of small Rq 
and 1/v [sij and spectra of stochastic fluctuations 
about each of them can be obtained separately 4J| . 
The theoretical analysis, however, does not allow us 
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Figure 5: (Color online) Simulation results for the (a), (d) dominant period of stochastic fluctuations, (b), 
(e) amplification, (c), (f) coherence and (g), (h) dominant period in the full spectrum. The seasonality is 
twice larger in the second row and panel (h) than in the first row and panel (g). The black stairstep graph 
bounds the gray region where all simulations went extinction within 500 years. Approximate parameter 
values for measles and rubella are depicted as dashed circles. Parameters: /i = 0.02 1/y, N = 10^, (a)-(c) 
and (g) e = 0.05, (d)-(f) and (h) e = 0.1. 
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to predict which of the attractors will be observed in 
simulations and what their relative contribution to 
the total stochastic dynamics is. Previous analysis 
of the stochastic dynamics of measles and pertussis 
showed that the only attractors seen in simulations 
of the seasonally forced SIR model (and other related 
models of infectious diseases) are annual and biennial 
cycles 0, S EHl • The stochastic switching was ob- 
served to happen exclusively between these attractors 
and only for measles. This result is, however, of lim- 
ited value to us because it is restricted to particular 
parameter choices, and so we cannot assume that the 
switching does not happen in the broader span of the 
parameter space. 

Spectra from simulations contain complete infor- 
mation about the frequency distribution of oscilla- 
tions and are thus helpful to identify switching be- 
tween attractors through the presence of unexpected 
peaks. Figure 5 shows simulation results for the dom- 
inant period of stochastic fluctuations, amplification 
and coherence for two values of seasonality e and 
= 10^. In addition to these quantities we com- 
pute the dominant period in the full spectrum which 
includes both stochastic and deterministic peaks. 

To examine the effect of seasonality on stochastic 
fluctuations Figure 5 (a)-(f) can be directly compared 
with Figure 3 (a)-(f) for the unforced model. As e in- 
creases the (gray) domain with frequent extinctions 
is extended and approaches the measles parameters. 
For most values of i?o and we have explored (the 
colored region) the dynamics of the stochastic model 
are associated with fluctuations about only few at- 
tractors. Firstly, the biennial cycle is found inside 
the region of increased coherence and amplification 
in Figure 5 (b)-(c) and (e)-(f) which includes measles 
and is absent in Figure 3 (b)-(c) and (c)-(f). The 
dominant period in the full spectrum of time series 
demonstrating such a behavior is at 2 years (Figure 
5 (h)). Secondly, stochastic switching between an- 
nual and triennial cycles is detected in a small region 
relevant for measles with low values of Rq, see an 
unexpected increase of amplification in Figure 5 (e) 
around i?o = 10. The spectra here have a domi- 
nant annual peak (Figure 5 (h)) and a subdominant 
triennial peak (Figure 5(d)). The amplification of os- 
cillations associated with the latter is however much 



higher than what we would expect to see for fluc- 
tuations around an annual cycle (see, e.g. Figure 5 
(b)). Thirdly, in the rest of the (colored) region which 
includes rubella, the spectra are similar to those of 
the unforced model. The dynamics of the stochas- 
tic model here corresponds to fluctuations about an 
annual cycle and we discuss it first. 

3.2.2 Implications for rubella and measles 
dynamics 

From Figure 5 (a) and (d) we see that for relatively 
small basic reproductive ratios, typical of rubella, the 
seasonality does not affect the dominant period of 
stochastic fluctuations which continues to be centered 
at about 5-6 years. The amplification (coherence) 
is only slightly increased (decreased) as e increases 
(Figure 5 (b)-(c) and (c)-(f)). The full spectra of 
rubella resemble that of Figure 3 with a sharp peak 
at 1 year and a broad multi-annual peak. 

For future discussion of the implications of vaccina- 
tion and decline or increase in birth rates, it is useful 
to investigate how the spectra of rubella change with 
Rq. Keeping the amplitude of seasonal forcing and 
infectious period of rubella fixed and increasing Rq , 
we expect the period of stochastic fiuctuations as well 
as their amplification to decrease. This is seen from 
Figure 5 and also illustrated in Figure 6 (a) where the 
full spectra for parameters close to rubella estimates 
are shown. The relative contribution of multi-annual 
and annual frequency components in model time se- 
ries can be read from the same figure. For small Rq 
the fiuctuations are large and the multi-annual peak 
is dominant but for larger Rq it becomes subdomi- 
nant and the annual peak is enhanced. The increase 
of e (or of population size as discussed before) re- 
sults in the enhancement of the deterministic peak 
too (compare Figure 5 (g) and (h)), but docs not 
change the dominant period of fluctuations signifi- 
cantly. 

For measles, there are major changes in the be- 
havior as both coherence and amplification increase 
drastically for e > 0, see the newly appeared oval- 
shaped regions near measles parameters in Figure 5 
(b)-(c) and (e)-(f). To demonstrate that this phe- 
nomenon indicates the appearance of a biennial cycle 
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(a) .... (b) 




Figure 6: (Color online) Full spectra from simulated time series, corresponding to the dynamics typical of (a) 
rubella and (b) measles. The color intensity increases (from blue to red) with increasing power. Parameters: 
H = 0.02 1/y, N = 10^ (a) e = 0.05, l/v = 18 d, (b) e = 0.1, l/i/ = 16 d. 



in simulations we show in Figure 6 (b) the full spec- 
tra for fixed infectious period, l/v = 16 d, and a 
range of Rq G [22,29]. These values of 1/z^ and Rq 
are slightly higher than the commonly used estimates 
for measles in England and Wales, e.g. 1/;^ = 13 d 
and i?o = 18 for London before vaccination. The 
same qualitative dynamics is observed for lower val- 
ues if the simulation length is shorter than 500 years 
(results not shown). For Rq ~ 22 the spectra are typ- 
ical of fluctuations about an annual limit cycle with 
a deterministic peak at 1 year and a broad stochas- 
tic peak near 2 years. If Rq is increased further the 
fluctuations around an annual cycle become macro- 
scopic (the stochastic peak at 2 years becomes much 
higher) smoothly turning into a biennial limit cycle. 
This transition corresponds to a period doubling bi- 
furcation in the deterministic model. A strong bien- 
nial behavior with a dominant peak at 2 years and a 
secondary harmonic at 1 year is observed for exam- 
ple, for Rq = 26. Finally, for even larger Rq we see 
a transition from biennial to an annual cycle again. 
The set of transitions seen in Figure 6 (b) is typical 
of measles and have been observed in related models 
of infections dynamics via analytical and numerical 
studies by other researches [H, O, [13, d, H, [HI]. 
The seasonality would act to change the picture in 



Figure 6 (b) in the following way. From the compar- 
ison of Figure 5 (g) and (h) the region of parameter 
space where such a behavior is seen is expected to get 
larger with increasing e, in particular for e > 0.1 the 
period doubling transition is induced for values of Rq 
much smaller than in Figure 6 (b). 

Previous analysis of measles data from England 
and Wales and the USA has shown that transitions 
in the dynamics due to an increase or decline of birth 
rates as well as the introduction of vaccination are 
associated with transition between annual and bien- 
nial limit cycles 17|. Using a simple mapping from 
changes in vaccination and birth rates to effective 
changes in Rq introduced in [l7| . our results confirm 
this view. For large communities with very high birth 
rates (high Rq) such as Liverpool before vaccination, 
US cities in the period after the Great Depression or 
cities in developing countries, we would expect to be 
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in the regime with an annual cycle [2^ [231, [17 
Other large cities with smaller birth rates such as 
London are in the regime with a biennial cycle [28 1. 
The corresponding spectrum with narrow and sharp 
peaks at 2 and 1 years has been the main reason of 
more regular and thus more predictable patterns of 
measles epidemics in large cities. The vaccination in- 
troduced in UK in 1968 lowered Rq and induced a 
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transition from the biennial to the annual cycle with 
large stochastic fluctuations. Our analysis thus offers 
an insight into the factors responsible for the shift 
from regular epidemics of measles before vaccination 
to less irregular in the vaccine era [43 1. 

The last finding deserving a further comment con- 
cerns the switching between an annual and triennial 
cycles found for moderate values of Rq, see Figure 5 
(e). This behavior may be responsible for the trien- 
nial cycles observed in Baltimore and other US cities 
during the Great Depression [s^ but more thorough 
analysis is needed to confirm this. 



4 Discussion 

In this paper, we have investigated the behavior of 
the stochastic seasonally forced SIR model based on 
spectra of long time series for a large range of basic 
reproductive ratios and infectious periods. For rela- 
tively low values of Rq relevant for rubella the model 
predicts spectra with a stochastic multi-annual peak 
at about 5-6 years and a deterministic annual peak. 
Both peaks are observed in the spectra of rubella data 
(Figure 1). The multi-annual peak stays largely un- 
changed under the introduction of seasonality (Fig- 
ures 3 and 5 (a), (d)) or population size (Figure 4) 
which explains its presence in time series from differ- 
ent locations. 

A visual inspection of simulated time series demon- 
strates intriguing behavior emerging from the interac- 
tion between stochasticity and a deterministic annual 
cycle. Figure 7 shows typical time series for the un- 
forced (red dashed line) and seasonally forced (blue 
solid line) cases. If e = the epidemic patterns rep- 
resent multi-annual coherent oscillations. As e is in- 
creased we find qualitatively different dynamics all of 
which correspond to spectra with a multi-annual and 
an annual peak. Figure 7 (a) shows an example of an- 
nual epidemics of alternating amplitudes modulated 
by an oscillation of a long period corresponding to the 
period of stochastic multi-annual fluctuations. This 
dynamic is qualitatively similar to the multi-annual 
regularity observed, for example, for rubella in Japan 
(Figure 1). We also find very large outbreaks followed 
by outbreaks of much lower amplitude as in Figure 7 



(b) . Such a behavior may be responsible for the spiky 
dynamics observed in, for example, Canada (Figure 
1). Note that the spikes in the data could also arise 
from spatial effects such as local extinction of the dis- 
ease followed by reintroduction from another region. 
However, as we do not possess more resolved data it 
is impossible to reach a final conclusion with regards 
to this issue. 

The patterns of rubella incidence in large popula- 
tions are in contrast with those of measles. For the 
latter the spectra are characterised by sharp and nar- 
row peaks at 1 and 2 years (as opposed to the broad 
multi-annual peak and a narrow peak at 1 year ob- 
served for rubella) and thus correspond to much more 
regular dynamics. The transitions in measles behav- 
ior due to vaccination or change in birth rates are 
associated with transitions between the annual and 
biennial limit cycles of the deterministic model. In fu- 
ture work, it would be interesting to study stochastic 
measles dynamics for higher levels of seasonal forcing 
that correspond to chaos in the deterministic model 
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Both measles and rubella are found to be close to 
the extinction boundary, and increasing the ampli- 
tude of seasonal forcing only extends the region of 
parameter space with frequent fade-outs. Figure 8 
illustrates typical stochastic trajectories from simu- 
lation in the susceptible-infected plane from which 
the spectra were computed. Interestingly, these pat- 
terns suggest that the mechanism accounting for high 
extinction probability arc different for rubella and 
measles. For rubella, extinction occurs as a conse- 
quence of large stochastic fluctuations about a small 
(and globally less stable) annual limit cycle, see Fig- 
ure 8 (a). For measles, extinctions are mainly due to 
the shape of a large (and globally more stable) bien- 
nial limit cycle, from which the system can be driven 
to extinction by even relatively small fluctuations, see 
Figure 8 (b). 

We can use the framework developed to predict the 
effect on persistence of an effective reduction in Eg 
by vaccination or declining birth rates for rubella and 
measles. For measles in the biennial regime, either an 
increase or a decrease of Rq can lead to fluctuations 
around an annual cycle (rather than a biennial cy- 
cle) that could result in lower extinction rates and 
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Figure 7: (Color online) Typical time series for rubella parameters. Parameters: N = 10®, /i = 0.02 1/y, 
i?o = 5, 1/zy = 18 d, e = (dashed red line), (a) e = 0.2 (solid blue line), (b) e = 0.3 (solid blue line). 




Figure 8: (Color online) Typical stochastic trajectories from simulation (blue solid lines) shown in the 
susceptible-infected plane. The red line is the deterministic annual cycle in (a) and the biennial cycle in (b). 
For the parameters used in this figure many simulations go extinct quite fast, for illustration purposes we 
have chosen the ones which lasted for 500 years. Parameters: = 10®, e = 0.05, fj, = 0.02 1/y, (a) Rq = 5, 
1/j/ = 18 d (close to rubella estimates), (b) Rq = 19, l/v = 12 A (close to measles estimates). 
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thus higher persistence. For rubella, a reduction in 
Rq will lead to larger and more coherent oscillations 
that would unambiguously result in higher extinc- 
tion probabilities, and thus lower persistence. Both 
these outcomes merit serious consideration in a pub- 
lic health context: vaccination against measles can 
make local elimination less likely [2l|; while vaccina- 
tion against rubella is likely to increase local extinc- 
tion, allowing the buildup of susceptible individuals 
in later age classes [s^, |33| , potentially leading to an 
increase in the burden of Congenital Rubella Syn- 
drome. 

These conclusions point to a need for theoretical 
developments towards uncovering the mechanisms of 
stochastic extinctions in small population based on 
the analysis of epidemic models (a thorough overview 
of studies in this area aiming at understanding the 
persistence of measles is given in the recent work by 



Conlan et. al. [l^)- In the mathematical framework 



we adopted in this paper, the approach to computa- 
tion of the distribution of extinction times in an un- 
forced stochastic epidemic model was proposed some 
Nevertheless, no analytical progress 



38, 3 



time ago 

can be made along the same lines for the seasonally 
forced model we use here. 

Our focus has been on measles and rubella; how- 
ever, the broad span of parameter space explored 
means that our results may shed light on the dy- 
namics of other diseases whose dynamics can be de- 
scribed by a simple SIR formalism with seasonal forc- 
ing; for example pertussis [i^ li^. Although the 
infection with pertussis does not confer permanent 
immunity the SIR model has been shown to capture 
the qualitative patterns to some extent [4^. Taking 
the pertussis parameters before vaccination that are 
well established in independent data sources (see e.g. 
0,11], 1/z/ = 22 d and i?o=17 for London and On- 
tario, Canada) from Figure 5 we see that for pertussis 
the dominant period of stochastic fluctuations is 2 to 
3 years. These periods are in the agreement with the 
documented interepidemic periods 
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We also 

find that the coherence of fluctuations is very low, 
which is compatible with the famously irregular dy- 
namics of pertussis. The decrease of i?o due to vac- 
cination would act to increase the dominant period 
and coherence which also agrees with the observation 



of the shift to more regular dynamics in the vaccine 
era @. 

To conclude, in this paper, we used a stochastic 
framework to explain the dynamics of rubella, 
particularly in comparison to measles. Our analysis 
revealed that whilst both rubella and measles are 
relatively close to their extinction boundary, the 
reasons for this are very different. Finally, our 
analysis showed that for rubella, reducing Rq by 
vaccinating, or a declining birth rate unambiguously 
result in higher extinction probabilities, whereas, for 
measles, outcomes can be more complicated; and 
both these facts have public health implications. 
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